# Comparison between Chop-Nod and Scan modes for Filament 5

# =============================================================================
# Package dependencies
# =============================================================================

from aplpy import FITSFigure
from astropy.wcs import WCS
from astropy.io import fits
from lmfit import Model
import matplotlib.pyplot as plt
import math
import numpy as np
from scipy import stats, optimize
from reproject import reproject_interp

# =============================================================================
# Magnetic field map on Stokes I total intensity
# =============================================================================
def BvImap(Iref, Pref, Region, FigSize, Ilevels, Iscale, IdI=10.0, PdP=3.0, Pmax=30.0,
           StepVec=3, ScaleVec=2, Colorbar='top', Linewidth=2.0, ScaleLength=None, PhysLength='1 pc'):
    # Plotting the INTENSITY POLARIZATION MAP
    BvI_plot = FITSFigure(Iref,figsize=(FigSize[0], FigSize[1]))
    BvI_plot.show_colorscale(cmap='Greys', vmin=Iscale[0], vmax=Iscale[1])
    BvI_plot.add_colorbar(pad=0.25)
    BvI_plot.colorbar.set_location(Colorbar)
    BvI_plot.colorbar.set_axis_label_text('Stokes $I$ Intensity (mJy arcsec$^{-2}$)')
    BvI_plot.add_beam(facecolor='white', edgecolor='red',
                        linewidth=2, pad=1, corner='bottom left') # HAWC+ Beam
    if ScaleLength != None:
        BvI_plot.add_scalebar(ScaleLength, PhysLength, corner='top right', frame=True) # Scalebar equivalent for 1 pc 
    # Adding total intensity contour
    BvI_plot.show_contour(data=Iref, colors='white', levels=Ilevels)
    # Centering
    BvI_plot.recenter(x=Region[0], y=Region[1], width=Region[2], height=Region[3])

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pref[8].data, header=Pref[0].header)
    pmap = pref.copy()
    # for the magnetic field angle B
    oref = fits.PrimaryHDU(data=Pref[11].data, header=Pref[0].header)
    omap = oref.copy()
    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pref[0].data/Pref[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pref[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pref[8].data/Pref[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pref[8].data > Pmax) # Polarization SNR threshold
    # Forcing the polarization vectors to share the same amplitude
    pmap.data[np.where(Pref[8].data > 0.0)] = 1.0
    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan

    # Plotting the polarization vectors
    BvI_plot.show_vectors(pmap, omap, scale=ScaleVec, step=StepVec, color='red', zorder=3, alpha=0.8, linewidth=Linewidth)

    return BvI_plot

# =============================================================================
# Polarization map on Stokes I total intensity
# =============================================================================
def PvImap(Iref, Pref, Region, FigSize, Ilevels, Iscale, IdI=10.0, PdP=3.0, Pmax=30.0,
           StepVec=3, ScaleVec=2, Colorbar='top', Linewidth=2.0):
    # Plotting the INTENSITY POLARIZATION MAP
    PvI_plot = FITSFigure(Iref,figsize=(FigSize[0], FigSize[1]))
    PvI_plot.show_colorscale(cmap='Greys', vmin=Iscale[0], vmax=Iscale[1])
    PvI_plot.add_colorbar(pad=0.25)
    PvI_plot.colorbar.set_location(Colorbar)
    PvI_plot.colorbar.set_axis_label_text('Stokes $I$ Intensity (mJy arcsec$^{-2}$)')
    PvI_plot.add_beam(facecolor='white', edgecolor='red',
                        linewidth=2, pad=1, corner='bottom left') # HAWC+ Beam
    # Adding the vector scale bar
    vectscale = ScaleVec * Iref.header['CDELT2'] # AplPy calculates vector scale with pixel units
    PvI_plot.add_scalebar(5 * vectscale, "P = 5%",corner='top right',frame=True) # Scalebar equivalent for 1 pc 
    # Adding total intensity contour
    PvI_plot.show_contour(data=Iref, colors='white', levels=Ilevels)
    # Centering
    PvI_plot.recenter(x=Region[0], y=Region[1], width=Region[2], height=Region[3])

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pref[8].data, header=Pref[0].header)
    pmap = pref.copy()
    # for the polarization angle O
    oref = fits.PrimaryHDU(data=Pref[10].data, header=Pref[0].header)
    omap = oref.copy()
    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pref[0].data/Pref[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pref[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pref[8].data/Pref[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pref[8].data > Pmax) # Polarization SNR threshold
    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan

    # Plotting the polarization vectors
    PvI_plot.show_vectors(pmap, omap, scale=ScaleVec, step=StepVec, color='red', zorder=3, alpha=0.8, linewidth=Linewidth)

    return PvI_plot

# =============================================================================
# Circular statistics
# =============================================================================
def circ_stat(array, base=1):
    # Takes a numpy array in degrees
    # Conversion from degrees to radians
    conv_rad = math.pi/180.0
    conv_deg = conv_rad**-1.0
    # Circular mean of the array
    circular_mean = conv_deg * stats.circmean(conv_rad*array, 
                                              low=base*-math.pi/2.0, 
                                              high=base*math.pi/2.0,
                                              nan_policy='omit')
    # Circular standard deviation of the array
    circular_stdev = conv_deg * stats.circstd(conv_rad*array, 
                                              low=base*-math.pi/2.0,  
                                              high=base*math.pi/2.0,
                                              nan_policy='omit')
    # Returning the values
    return circular_mean, circular_stdev


# =============================================================================
# Function to create an anngle comparison plot between Chop-Nod and Scan mode data
# =============================================================================
def CompareAngles(RefFit, ProjFit, FigSize=[5.12,3.84], IdI=10, PdP=3, Pmax=30):
    # Creating working numpy arrays for the angles
    RefArrayO= np.copy(RefFit[10].data)
    ProjArrayO = np.copy(ProjFit[10].data)
    RefArraydO= np.copy(RefFit[12].data)
    ProjArraydO = np.copy(ProjFit[12].data)
    # Creating the lists of angles for each map
    RefCatO = []
    ProjCatO = []
    RefCatdO = []
    ProjCatdO = []

    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    # Reference cube
    imask_01 = np.where(RefFit[0].data/RefFit[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(RefFit[0].data < 0) # Total intensity positive threshold
    pmask01 = np.where(RefFit[8].data/RefFit[9].data < PdP) # Polarization SNR threshold
    pmaxmask01 = np.where(RefFit[8].data > Pmax) # Polarization SNR threshold
    # Projected cube
    imask_03 = np.where(ProjFit[0].data/ProjFit[1].data < IdI) # Total intensity SNR threshold
    imask_04 = np.where(ProjFit[0].data < 0) # Total intensity positive threshold
    pmask02 = np.where(ProjFit[8].data/ProjFit[9].data < PdP) # Polarization SNR threshold
    pmaxmask02 = np.where(ProjFit[8].data > Pmax) # Polarization SNR threshold
    # Applying all the SNR masks to the reference array
    RefArrayO[imask_01] = np.nan
    RefArrayO[imask_02] = np.nan
    RefArrayO[pmask01] = np.nan
    RefArrayO[pmaxmask01] = np.nan
    RefArrayO[imask_03] = np.nan
    RefArrayO[imask_04] = np.nan
    RefArrayO[pmask02] = np.nan
    RefArrayO[pmaxmask02] = np.nan

    # Measuring the size of the array
    shapeI = RefArrayO.shape[0]
    shapeJ = RefArrayO.shape[1]
    # Creating the catalogs
    for i in range (0,shapeI):
        for j in range (0,shapeJ):
            if RefArrayO[i,j] > -360.0: # Iterate only on pixels that a,.re not NaN
                RefCatO.append(RefArrayO[i,j])
                ProjCatO.append(ProjArrayO[i,j])
                RefCatdO.append(RefArraydO[i,j])
                ProjCatdO.append(ProjArraydO[i,j])
    # Catalog of differences
    DiffCat = np.array(ProjCatO) - np.array(RefCatO) # Not robust, this only works here because no differences are greater than 90 degree
    #DiffCatAbs = np.abs(np.array(ProjCatO) - np.array(RefCatO))
    # Statistical properties
    MeanO,StdevO = circ_stat(DiffCat)
    print('Differences between catalogs')
    print('Circular mean: '+str(MeanO))
    print('Circular standard deviation: '+str(StdevO))
    print('Maximum: ' + str(np.max(np.array(DiffCat))))
    print('Minimum: ' + str(np.min(np.array(DiffCat))))
    # Median uncertainties
    print('Mean angle uncertainty (Chop-Nod, Scan): (' + str(np.mean(np.array(RefCatdO))) + ',' + str(np.mean(np.array(ProjCatdO))) +')')
    print('Median angle uncertainty (Chop-Nod, Scan): (' + str(np.median(np.array(RefCatdO))) + ',' + str(np.median(np.array(ProjCatdO))) +')')
    print('Standard deviation angle uncertainty (Chop-Nod, Scan): (' + str(np.std(np.array(RefCatdO))) + ',' + str(np.std(np.array(ProjCatdO))) +')')
    print('Max angle uncertainty (Chop-Nod, Scan): (' + str(np.max(np.array(RefCatdO))) + ',' + str(np.max(np.array(ProjCatdO))) +')')
    print('Min angle uncertainty (Chop-Nod, Scan): (' + str(np.min(np.array(RefCatdO))) + ',' + str(np.min(np.array(ProjCatdO))) +')')
    # MeanO_Abs,StdevO_Abs = circ_stat(DiffCatAbs)
    # print('Absolute differences between catalogs')
    # print('Circular mean: '+str(MeanO_Abs))
    # print('Circular standard deviation: '+str(StdevO_Abs))

    # Creating the comparison figure
    CompA_plot = plt.figure(figsize=(FigSize[0],FigSize[1]))
    plt.plot(RefCatO, ProjCatO, '+', color='black', markersize=4)
    # Zero lines
    plt.axvline(x=0.0, color='black', linestyle='--', linewidth=1)
    plt.axhline(y=0.0, color='black', linestyle='--', linewidth=1)
    # Creating zero difference line
    RefX = np.arange(-90.0, 90.0, 0.1)
    RefY = np.arange(-90.0, 90.0, 0.1)
    plt.plot(RefX, RefY, 'black')
    # Creating difference lines
    Ref45Xa = np.arange(-90.0, 45.0, 0.1)
    Ref45Ya = np.arange(-45.0, 90.0, 0.1)
    plt.plot(Ref45Xa, Ref45Ya, 'black', linestyle=':', linewidth=1)
    Ref45Yb = np.arange(-90.0, 45.0, 0.1)
    Ref45Xb = np.arange(-45.0, 90.0, 0.1)
    plt.plot(Ref45Xb, Ref45Yb, 'black', linestyle=':', linewidth=1)
    Ref45Xc = np.arange(-90.0, -45.0, 0.1)
    Ref45Yc = np.arange(45.0, 90.0, 0.1)
    plt.plot(Ref45Xc, Ref45Yc, 'black', linestyle=':', linewidth=1)
    Ref45Yd = np.arange(-90.0, -45.0, 0.1)
    Ref45Xd = np.arange(45.0, 90.0, 0.1)
    plt.plot(Ref45Xd, Ref45Yd, 'black', linestyle=':', linewidth=1)
    Ref90Xa = np.arange(-90.0, 0.0, 0.1)
    Ref90Ya = np.arange(0.0, 90.0, 0.1)
    plt.plot(Ref90Xa, Ref90Ya, 'black', linestyle='-.', linewidth=1)
    Ref90Yb = np.arange(-90.0, 0.0, 0.1)
    Ref90Xb = np.arange(0.0, 90.0, 0.1)
    plt.plot(Ref90Xb, Ref90Yb, 'black', linestyle='-.', linewidth=1)
    # Adding the labels to the difference lines
    plt.text(-75, -85, r'$\Delta \theta = 0^{\circ}$')
    plt.text(-30, -85, r'$\Delta \theta = 45^{\circ}$')
    plt.text(15, -85, r'$\Delta \theta = 90^{\circ}$')
    plt.text(60, -85, r'$\Delta \theta = 45^{\circ}$')
    plt.text(-80, 80, r'$\Delta \theta = 45^{\circ}$')
    plt.text(-35, 80, r'$\Delta \theta = 90^{\circ}$')
    plt.text(10, 80, r'$\Delta \theta = 45^{\circ}$')
    plt.text(55, 80, r'$\Delta \theta = 0^{\circ}$')

    # Setting axes and labels
    plt.ylabel(r'Scan Mode Polarization Angle (degree)')
    plt.xlabel(r'Chop-Nod Polarization Angle (degree)')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side   
    plt.xlim(-90.0,90.0) # X axis range
    plt.ylim(-90.0,90.0) # Y axis range
    plt.tight_layout() # Using all available space in the plot window

    return CompA_plot

# =============================================================================
# Function to create an angle comparison map between Chop-Nod and Scan mode data
# =============================================================================
def CompareAnglesMap(Iref, RefFit, ProjFit, Region, FigSize, Ilevels, IdI=10, PdP=3, Pmax=30, Colorbar='top'):
    # Creating working numpy arrays for the angles
    RefArrayO= np.copy(RefFit[10].data)
    ProjArrayO = np.copy(ProjFit[10].data)

    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    # Reference cube
    imask_01 = np.where(RefFit[0].data/RefFit[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(RefFit[0].data < 0) # Total intensity positive threshold
    pmask01 = np.where(RefFit[8].data/RefFit[9].data < PdP) # Polarization SNR threshold
    pmaxmask01 = np.where(RefFit[8].data > Pmax) # Polarization SNR threshold
    # Projected cube
    imask_03 = np.where(ProjFit[0].data/ProjFit[1].data < IdI) # Total intensity SNR threshold
    imask_04 = np.where(ProjFit[0].data < 0) # Total intensity positive threshold
    pmask02 = np.where(ProjFit[8].data/ProjFit[9].data < PdP) # Polarization SNR threshold
    pmaxmask02 = np.where(ProjFit[8].data > Pmax) # Polarization SNR threshold
    # Applying all the SNR masks to the reference array
    RefArrayO[imask_01] = np.nan
    RefArrayO[imask_02] = np.nan
    RefArrayO[pmask01] = np.nan
    RefArrayO[pmaxmask01] = np.nan
    RefArrayO[imask_03] = np.nan
    RefArrayO[imask_04] = np.nan
    RefArrayO[pmask02] = np.nan
    RefArrayO[pmaxmask02] = np.nan

    # Measuring the size of the array
    shapeI = RefArrayO.shape[0]
    shapeJ = RefArrayO.shape[1]
    # Creating the difference map
    DiffMapO =np.empty((shapeI,shapeJ))
    DiffMapO[:] = np.nan
    for i in range (0,shapeI):
        for j in range (0,shapeJ):
            if RefArrayO[i,j] > -360.0: # Iterate only on pixels that are not NaN
                Diff_Value = RefArrayO[i,j] - ProjArrayO[i,j]
                # Differences should never be more than -90 or 90 degrees
                if Diff_Value > 90.0:
                    Diff_Value = Diff_Value - 180.0
                if Diff_Value < -90.0:
                    Diff_Value = Diff_Value + 180.0
                # Adding value to the difference map
                DiffMapO[i,j] = Diff_Value
    # Creating fits container
    DiffRef = fits.PrimaryHDU(data=DiffMapO, header=RefFit[0].header)
    # Plotting the resulting map
    CompA_map = FITSFigure(DiffRef,figsize=(FigSize[0], FigSize[1]))
    CompA_map.show_colorscale(cmap='rainbow', vmin=-50.0, vmax=50.0)
    CompA_map.add_colorbar(pad=0.25)
    CompA_map.colorbar.set_location(Colorbar)
    CompA_map.colorbar.set_axis_label_text(r'Polarization Angle Difference $\theta_{ch} - \theta_{sc}$ (degree)')
    CompA_map.add_beam(facecolor='white', edgecolor='red',
                        linewidth=2, pad=1, corner='bottom left') # HAWC+ Beam
    # Adding total intensity contour
    CompA_map.show_contour(data=Iref, colors='black', levels=Ilevels)
    # Centering
    CompA_map.recenter(x=Region[0], y=Region[1], width=Region[2], height=Region[3])

    return CompA_map

# =============================================================================
# Function to create an intensity comparison plot between Chop-Nod and Scan mode data
# =============================================================================
def CompareIntensity(RefFit, ProjFit, Iscale=[0,500], FigSize=[5.12,3.84]):
    # Conversion factor from Jy/pixel to mJy/arcsec^2
    fluxconv = 1000.0*(3600.0*RefFit[0].header['CDELT2'])**-2.0
    # Creating the lists of Stokes I for each map
    RefCat = []
    ProjCat = []
    
    # Creating working numpy arrays
    RefArray = np.copy(RefFit[0].data)  
    ProjArray = np.copy(ProjFit[0].data)
    # Masking all negative values for the total intensity
    imask_01 = np.where(RefFit[0].data < 0)
    imask_02 = np.where(ProjFit[0].data < 0)
    # Masking the data in the reference array
    RefArray[imask_01] = np.nan
    RefArray[imask_02] = np.nan
    
    # Measuring the size of the array
    shapeI = RefArray.shape[0]
    shapeJ = RefArray.shape[1]
    # Creating the catalogs
    for i in range (0,shapeI):
        for j in range (0,shapeJ):
            if RefArray[i,j] > 0.0: # Iterate only on pixels that are not NaN
                RefCat.append(RefArray[i,j]* fluxconv)
                ProjCat.append(ProjArray[i,j]* fluxconv) 
    
    # Creating reference line
    linemin = Iscale[0]
    linemax = Iscale[1]
    RefX = np.arange(linemin, linemax, 0.1)
    RefY = np.arange(linemin, linemax, 0.1)
    
    # Fitting a linear relation to IvI
    def IvIfunc(x, a, b):
        return a*x+b
    IvIpara, IvIcov = optimize.curve_fit(IvIfunc, RefCat, ProjCat)
    # Printing the results
    print('====')
    print('Fit parameters for y = a*x+b')
    print('Coefficient: '+str(IvIpara[0]))
    print('Constant: '+str(IvIpara[1]))
    print('====')
    # Creating x values
    IvIfuncX = np.arange(linemin, linemax, 0.1)
    # Calculating the y values
    IvIfuncY = IvIfunc(IvIfuncX,IvIpara[0],IvIpara[1])

    # Constructing the plot
    IvI_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Chop-Nod Stokes I Flux (mJy arcsec$^{-2}$)')
    plt.ylabel('Scan Stokes I Flux (mJy arcsec$^{-2}$)')
    plt.plot(RefCat, ProjCat, '.', color='black', markersize=1, mfc='none') # Creating the data points
    plt.plot(RefX, RefY, 'orange') # Creating the 1:1 reference line
    plt.plot(IvIfuncX, IvIfuncY, 'red') # Plotting the fitted linear relation
    # Modifying the axes
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    axes.set_xlim(linemin, linemax)
    axes.set_ylim(linemin, linemax)
    plt.tight_layout() # Using all available space in the plot window
    return IvI_plot

# =============================================================================
# Function to create an intensity comparison map between Chop-Nod and Scan mode data
# =============================================================================
def CompareIntensityMap(Iref, RefFit, ProjFit, Region, FigSize, Ilevels, IdI=10, PdP=3, Pmax=30, Colorbar='top'):
    # Creating working numpy arrays
    RefArray = np.copy(RefFit[0].data)  
    ProjArray = np.copy(ProjFit[0].data)
    # Masking all negative values for the total intensity
    imask_01 = np.where(RefFit[0].data < 0)
    imask_02 = np.where(ProjFit[0].data < 0)
    # Masking the data in the reference array
    RefArray[imask_01] = np.nan
    RefArray[imask_02] = np.nan

    # Measuring the size of the array
    shapeI = RefArray.shape[0]
    shapeJ = RefArray.shape[1]
    # Creating the difference map
    RatioMap=np.empty((shapeI,shapeJ))
    RatioMap[:] = np.nan
    # Creating the catalogs
    for i in range (0,shapeI):
        for j in range (0,shapeJ):
            if RefArray[i,j] > 0.0: # Iterate only on pixels that are not NaN
                RatioMap[i,j] = ProjArray[i,j]/RefArray[i,j]

    # Creating fits container
    RatioRef = fits.PrimaryHDU(data=RatioMap, header=RefFit[0].header)
    # Plotting the resulting map
    CompI_map = FITSFigure(RatioRef,figsize=(FigSize[0], FigSize[1]))
    CompI_map.show_colorscale(cmap='rainbow', vmin=0.0, vmax=2.0)
    CompI_map.add_colorbar(pad=0.25)
    CompI_map.colorbar.set_location(Colorbar)
    CompI_map.colorbar.set_axis_label_text(r'Intensity Ratio $I_{sc} I^{-1}_{ch}$')
    CompI_map.add_beam(facecolor='white', edgecolor='red',
                        linewidth=2, pad=1, corner='bottom left') # HAWC+ Beam
    # Adding total intensity contour
    CompI_map.show_contour(data=Iref, colors='black', levels=Ilevels)
    # Centering
    CompI_map.recenter(x=Region[0], y=Region[1], width=Region[2], height=Region[3])

    return CompI_map

# =============================================================================
# Function to create an Histogram and fit a Gaussian function
# =============================================================================
def Histogram(CatdO, binsize=5.0, showfit='no', FigSize=[5.12,3.0], Source=None, LabelSide='left', FilO=None, PlanckO=None):
    # Size of the catalog
    CatSize = CatdO.shape[0]
    
    # Method to create an histogram from a vector catalog
    print('Calculating mean and standard deviation of the catalog')
    print()
    
    # Conversion from degrees to radians
    conv_rad = math.pi/180.0
    conv_deg = conv_rad**-1.0
    
    # Calculating the regular mean
    norm_mean = np.mean(CatdO)
    norm_std = np.std(CatdO)
    
    print('Mean: ' + str(norm_mean))
    print('Standard deviation: ' + str(norm_std))
    print()
    
    # Calculating the circular mean and circular standard deviation
    # assuming boundaries of -90 to 90 degrees
    circ_mean = conv_deg * stats.circmean(conv_rad*CatdO, 
                                            low=-math.pi/2.0, 
                                            high=math.pi/2.0)
    circ_std = conv_deg * stats.circstd(conv_rad*CatdO, 
                                        low=-math.pi/2.0,  
                                        high=math.pi/2.0)
    
    print('Circular mean: ' + str(circ_mean))
    print('Circular standard deviation: ' + str(circ_std))
    print()
    
    # Calculating the number of bins by rounding to lowest integer
    number_bins = math.floor(180.0 / (1.0 * binsize))
    eff_binsize = 180.0/number_bins
    
    print('Number of bins: ' + str(number_bins))
    print('Bin size used: ' + str(eff_binsize))
    print()
    
    # Creating the histogram arrays
    histo_bins = np.zeros(number_bins)
    histo_values = np.zeros(number_bins)
    # Identifying the values of the bins
    for i in range (0,number_bins):
        histo_bins[i] = (i + 0.5) * eff_binsize - 90.0
    # Counting the number of vectors per bin - BRUTE FORCE VERSION
    # CONSIDER USING NUMPY HISTOGRAM FUNCTION INSTEAD
    for i in range (0,number_bins):
        # Limits of the bin
        min_limit = histo_bins[i] - 0.5 * eff_binsize
        max_limit = histo_bins[i] + 0.5 * eff_binsize
        # Checking every vector if they should be in the bin
        for j in range (0,CatSize):
            if (CatdO[j] >= min_limit) and (CatdO[j] < max_limit):
                histo_values[i] = histo_values[i] + 1
    
    # Find array value closest to circular mean    
    abs_array = np.abs(histo_bins - circ_mean)
    smallest_diff = abs_array.argmin()        
    # Calculating shift needed to center the histogram on circular mean
    shift = math.ceil(number_bins/2.0) - smallest_diff
    # Rolling the histogram
    print('Using the circular mean to center the histogram')
    print('Rolling histogram by '+str(shift)+' elements')
    print()
    histo_bins_rolled = np.roll(histo_bins, shift)
    histo_values_rolled = np.roll(histo_values, shift)
    
    # Fixing values of rolled bins so array stays sorted
    if shift > 0:
        for i in range(0, shift):
            histo_bins_rolled[i] = histo_bins_rolled[i] - 180.0
    if shift < 0:
        for i in range(number_bins+shift, number_bins):
            histo_bins_rolled[i] = histo_bins_rolled[i] + 180.0
    
    # Attempting Gaussian fit to the data
    # Defining the Gaussian function, stolen from the internet
    # https://lmfit.github.io/lmfit-py/model.html
    def gaussian(x, amp, cen, wid):
        y = amp*np.exp(-(x-cen)**2 / (2*wid**2))
        return y
    
    print('Fitting a Gaussian profile to the histogram')
    print()
    
    # Calling the lmfit package to model the 'gaussian' function                        
    gaussian_model = Model(gaussian)
    gaussian_fit = gaussian_model.fit(histo_values_rolled, 
                                        x=histo_bins_rolled, 
                                        amp=np.mean(histo_values_rolled), 
                                        cen=circ_mean, wid=circ_std)
    # Creating the array containing the parameters
    fit_params = np.empty([3,2])
    # Amplitude
    fit_params[0,0] = gaussian_fit.params['amp'].value
    fit_params[0,1] = gaussian_fit.params['amp'].stderr
    # Mean
    fit_params[1,0] = gaussian_fit.params['cen'].value
    fit_params[1,1] = gaussian_fit.params['cen'].stderr
    # Standard deviation
    fit_params[2,0] = gaussian_fit.params['wid'].value
    fit_params[2,1] = gaussian_fit.params['wid'].stderr

    print('Goodness of fit')
    print('Number of iterations: '+str(gaussian_fit.nfev))
    print('Reduced chi-squared: '+str(gaussian_fit.redchi)) 
    print()
    
    print('Fitted parameters')
    print('Amplitude:          '+str(gaussian_fit.params['amp'].value))
    print('                  ± '+str(gaussian_fit.params['amp'].stderr))
    print()
    print('Mean:               '+str(gaussian_fit.params['cen'].value))
    print('                  ± '+str(gaussian_fit.params['cen'].stderr))
    print()
    print('Standard deviation: '+str(gaussian_fit.params['wid'].value))
    print('                  ± '+str(gaussian_fit.params['wid'].stderr))
    print()
    
    # Plotting the histogram
    histo_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.bar(histo_bins_rolled, histo_values_rolled, 
                    width=eff_binsize, fill=False) # Creating the bar plot
    plt.xlabel(r'Polarization Angle $\theta$ (Degree)')
    plt.ylabel('Number')
    plt.tight_layout() # Using all available space in the plot window
    if Source != None:
        # Creating the reference box
        props = dict(boxstyle='square', facecolor='white')
        # Including the information
        textstr = '\n'.join((
            Source,
            r'$\overline{\theta}=%.1f^{\circ}$' % circ_mean,
            r'$\sigma_{\theta}=%.1f^{\circ}$'% circ_std))
        if FilO != None:
            textstr = '\n'.join((
            Source,
            r'$\overline{\theta}=%.1f^{\circ}$' % circ_mean,
            r'$\sigma_{\theta}=%.1f^{\circ}$'% circ_std,
            r'$\phi_{F}=%.1f^{\circ}$'% FilO))
        if PlanckO != None:
            textstr = '\n'.join((
            Source,
            r'$\overline{\theta}=%.1f^{\circ}$' % circ_mean,
            r'$\sigma_{\theta}=%.1f^{\circ}$'% circ_std,
            r'$\overline{\theta}_{P}=%.1f^{\circ}$'% PlanckO))
        if FilO != None and PlanckO !=None:
            textstr = '\n'.join((
            Source,
            r'$\overline{\theta}=%.1f^{\circ}$' % circ_mean,
            r'$\sigma_{\theta}=%.1f^{\circ}$'% circ_std,
            r'$\overline{\theta}_{P}=%.1f^{\circ}$'% PlanckO,
            r'$\phi_{F}=%.1f^{\circ}$'% FilO))
        if LabelSide == 'left':
            histo_plot.text(0.17, 0.75, textstr, verticalalignment='center', horizontalalignment='left', bbox=props, fontweight='bold')
        if LabelSide == 'right':
            histo_plot.text(0.94, 0.75, textstr, verticalalignment='center', horizontalalignment='right', bbox=props, fontweight='bold')
    # Setting range
    xmin = histo_bins_rolled[0]-eff_binsize/2.0
    xmax = histo_bins_rolled[number_bins-1]+eff_binsize/2.0
    plt.xlim(xmin, xmax)
    # Modifying the axes directly
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    # Plotting the best fit Gaussian
    # Creating x values
    gaussian_bins = np.arange(xmin, xmax, 0.1)
    # Calculating the y values
    gaussian_values = gaussian(gaussian_bins, fit_params[0,0], 
                                fit_params[1,0], fit_params[2,0])
    # Plotting the Gaussian fit
    if showfit == 'yes':
        plt.plot(gaussian_bins, gaussian_values, 'k')

    # Adding reference lines
    plt.axvline(x=-180.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=-90.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=0.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=90.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=180.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=circ_mean, color='red', linestyle='-', linewidth=2)
    if FilO != None: 
        plt.axvline(x=FilO, color='blue', linestyle='-.', linewidth=2)
    if PlanckO != None: 
        plt.axvline(x=PlanckO, color='green', linestyle='-.', linewidth=2)
    
    # Returning figure and fits parameters
    return histo_plot, fit_params

# =============================================================================
# Creating the polarization vector catalog
# =============================================================================
def MakeCats(Iref, Pref, CatMask, IdI=10.0, PdP=3.0, Pmax=30.0):
    Idata = Iref.copy() # Total intensity map converted to the correct units
    Pdata = Pref.copy() # Full Polarization data cube
    Maskdata = CatMask.copy() # Mask for the catalog

    # Creating lists
    CatI = [] # Stokes I
    CatdI =[] # Uncertainty dI for Stokes I
    CatQ =[] # Stokes Q
    CatdQ =[] # Uncertainty dQ for Stokes Q
    CatU =[] # Stokes U
    CatdU =[] # Uncertainty dU for Stokes U
    CatP = [] # Polarization fraction P (debiased)
    CatdP =[] # Uncertainty dP for polarization fraction P
    CatO = [] # Polarization angle O
    CatB = [] # Rotated polarization angle B (+90 degrees)
    CatdO =[] # Uncertainty dO for polarization angles
    CatPI = [] # Polarized intensity PI (debiased)
    CatdPI =[] # Uncertainty dPI for polarized intensity PI

    fluxconv = 1000.0*(3600.0*Idata.header['CDELT2'])**-2.0 # Conversion factor from Jy/pixel

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pdata[8].data, header=Pdata[0].header)
    pmap = pref.copy()
    # for the magnetic field angle B
    oref = fits.PrimaryHDU(data=Pdata[10].data, header=Pdata[0].header)
    omap = oref.copy()
    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pdata[0].data/Pdata[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pdata[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pdata[8].data/Pdata[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pdata[8].data > Pmax) # Polarization SNR threshold
    cmask = np.where(Maskdata.data != 1.0) # Catalog mask
    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan
    pmap.data[cmask] = np.nan

    # Measuring the size of the new array
    shapeX = Pdata[0].data.shape[0]
    shapeY = Pdata[0].data.shape[1]

    for i in range (0,shapeX):
        for j in range (0,shapeY):
            if pmap.data[i,j] > 0.0: # Iterate only on pixels where the selection criteria are respected
                CatI.append(Idata.data[i,j])
                CatdI.append(Pdata[1].data[i,j]* fluxconv)  # Conversion factor from Jy/pixel
                CatQ.append(Pdata[2].data[i,j]* fluxconv)
                CatdQ.append(Pdata[3].data[i,j]* fluxconv)
                CatU.append(Pdata[4].data[i,j]* fluxconv)
                CatdU.append(Pdata[5].data[i,j]* fluxconv)
                CatP.append(pmap.data[i,j])
                CatdP.append(Pdata[9].data[i,j])
                CatO.append(Pdata[10].data[i,j])
                CatB.append(Pdata[11].data[i,j])
                CatdO.append(Pdata[12].data[i,j])
                CatPI.append(Pdata[15].data[i,j]* fluxconv)
                CatdPI.append(Pdata[14].data[i,j]* fluxconv)

    # Transforming lists into numpy arrays
    CatI = np.array(CatI)
    CatdI = np.array(CatdI)
    CatQ = np.array(CatQ)
    CatdQ = np.array(CatdQ)
    CatU = np.array(CatU)
    CatdU = np.array(CatdU)
    CatP = np.array(CatP)
    CatdP = np.array(CatdP)
    CatO = np.array(CatO)
    CatB = np.array(CatB)
    CatdO = np.array(CatdO)
    CatPI = np.array(CatPI)
    CatdPI = np.array(CatdPI)

    # Returning all the catalogs 
    return CatI, CatdI, CatQ, CatdQ, CatU, CatdU, CatP, CatdP, CatO, CatB,CatdO, CatPI, CatdPI

# =============================================================================
# Creating a mask for a region of the map
# =============================================================================
def MaskRegion(FitsFile, CatRegion, angle=0):
       # Reading region information from CatRegion
       center_RA = CatRegion[0]
       center_DEC = CatRegion[1]
       width = CatRegion[2]
       height = CatRegion[3]
       
       # Creating the empty array for the final mask
       maskref = fits.PrimaryHDU(data=FitsFile[0].data, header=FitsFile[0].header)
       mask = maskref.copy()
       # Size of the map
       MaskSize = np.shape(FitsFile[0].data) 

       # Creating WCS object
       MaskWCS = WCS(FitsFile[0].header)          
       # Identifying pixels to be included in the mask
       for i in range (0,MaskSize[0]): # Declination
              for j in range (0,MaskSize[1]): # Right Ascension
                     # Determining the current ICRS position in the array
                     MaskRA, MaskDec = (MaskWCS.array_index_to_world_values(i, j))
                     # Distance to central pixel
                     delta_ra = (MaskRA - 
                                   center_RA) * math.cos(
                                   MaskDec*math.pi/180.0)
                     delta_dec = MaskDec - center_DEC

                     # Coordinate conversion
                     x1 = delta_ra * -1.0
                     y1 = delta_dec
                     ang = angle*math.pi/180.0
                     # Measuring angular distances
                     delta_x2 = abs(x1*math.cos(ang)+y1*math.sin(ang))
                     delta_y2 = abs(-x1*math.sin(ang)+y1*math.cos(ang))
                     # Is the pixel in the mask? If yes, set to 1.0, if no, set to NaN
                     if delta_x2 < 0.5*width and delta_y2 < 0.5*height:
                            mask.data[i,j] = 1.0
                     else:
                            mask.data[i,j] = np.nan
       return mask

# =============================================================================
# Polarization fraction as a function of total intensity
# =============================================================================
def PvI(CatP, CatdP, CatI, Iscale=None, Pscale=None, showfit='false', weighted='true', FigSize=[5.12,3.84], Source=None):  
    PvI_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Stokes Total Intensity $I$ (mJy arcsec$^{-2}$)')
    plt.ylabel("Polarization Fraction $P$ (%)")
    
    plt.plot(CatI, CatP, '.', color='black', markersize=1, mfc='none')
    plt.xscale('log')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    plt.tight_layout() # Using all available space in the plot window

    # Checking the automated limits for I and P
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    # Setting mininum and maximum values of the I plot
    if Iscale==None:
        Iscale = np.empty(2)
        Iscale[0] = xmin
        Iscale[1] = xmax
    if Pscale==None:
        Pscale = np.empty(2)
        Pscale[0] = ymin
        Pscale[1] = ymax
    # Setting the minimum and maximum values
    plt.xlim(Iscale[0],Iscale[1])
    plt.ylim(Pscale[0],Pscale[1])


    if showfit=='true':
        # Create the array with ln(I), ln(P), and d(ln(P))
        # The numpy package uses log(x) as the natural logarithm ln
        length = CatP.size
        PvIarray = np.empty([length, 3])
        PvIarray[:,0]=np.log(CatI)
        PvIarray[:,1]=np.log(CatP)
        PvIarray[:,2]=CatdP/CatP
        # Defining linear function to be fitted by curve_fit
        def lnPvlnI(x, a, b):
            return a*x + b
        # Fitting power law to data
        if weighted=='true':
            PvIpow, PvIcov = optimize.curve_fit(lnPvlnI, PvIarray[:,0], PvIarray[:,1], sigma=PvIarray[:,2])
        else:
            PvIpow, PvIcov = optimize.curve_fit(lnPvlnI, PvIarray[:,0], PvIarray[:,1])
        # Uncertainties on the fit
        PvIerr = np.sqrt(np.diag(PvIcov))
        # Converting the coefficient back to its non-logarithmic value
        PvIcoeff = math.exp(PvIpow[1])
        PvIdcoeff = math.exp(PvIpow[1]) * PvIerr[1]
        # Printing the results of the fit
        print('Power law index: '+str(PvIpow[0])+' ± '+str(PvIerr[0]))
        print('Coefficient: '+str(PvIcoeff)+' ± '+str(PvIdcoeff))
        # Calculating the reduced Chi-Squared value
        if weighted=='false':
            FakeData = lnPvlnI(PvIarray[:,0],PvIpow[0],PvIpow[1])
            Chisq = np.sum((PvIarray[:,0]-FakeData)**2)
            print('Chi-Squared: ' + str(Chisq))
            print('Number of elements: ' + str(PvIarray[:,0].size))
            RedChisq = Chisq/(PvIarray[:,0].size-2) # Reduced Chi-Squared is number of elements minus number of fitted parameters
            print('Reduced Chi-Squared: ' +str(RedChisq))

        if Source != None:
            PvI_plot.text(0.92, 0.85, Source, verticalalignment='bottom', horizontalalignment='right', fontweight='bold')

    # Adding power law fit
    if showfit=='true':
        def PvIfunc(x, A, B):
            return A*x**B
            
        # Creating x values
        PvIfuncX = np.arange(Iscale[0], Iscale[1], 0.1)
        # Calculating the y values
        PvIfuncY = PvIfunc(PvIfuncX,PvIcoeff,PvIpow[0])
        plt.plot(PvIfuncX, PvIfuncY, 'red')
        # Rounding up the fit error
        RoundErrorPvI = round_decimals_up(PvIerr[0],2)
        RoundErrorCoeff = round_decimals_up(PvIdcoeff,1)

        # Adding the fit results as a label
        Coefficients = '\n'.join((r'$\alpha = {0:.2f} \pm {1:.2f}$'.format(PvIpow[0],RoundErrorPvI),
                                 r'$A = {0:.1f} \pm {1:.1f}$'.format(PvIcoeff,RoundErrorCoeff)))

        PvI_plot.text(0.92, 0.75, Coefficients, verticalalignment='bottom', horizontalalignment='right')

    return PvI_plot

# =============================================================================
# Function to reproject and re-calculate a polarization cube
# =============================================================================
def ReprojectCube(RefFit, TempFit):
    # Reprojecting the Tempfit data to the correct pixel scale
    ProjI, footprint = reproject_interp(TempFit[0], RefFit[0].header) # Stokes I
    ProjI_fit = fits.ImageHDU(data=ProjI, header=RefFit[0].header) 
    ProjdI, footprint = reproject_interp(TempFit[1], RefFit[1].header) # Stokes dI
    ProjdI_fit = fits.ImageHDU(data=ProjdI, header=RefFit[1].header) 
    ProjQ, footprint = reproject_interp(TempFit[2], RefFit[2].header) # Stokes Q
    ProjQ_fit = fits.ImageHDU(data=ProjQ, header=RefFit[2].header)
    ProjdQ, footprint = reproject_interp(TempFit[3], RefFit[3].header) # Stokes dQ
    ProjdQ_fit = fits.ImageHDU(data=ProjdQ, header=RefFit[3].header)
    ProjU, footprint = reproject_interp(TempFit[4], RefFit[4].header) # Stokes U
    ProjU_fit = fits.ImageHDU(data=ProjU, header=RefFit[4].header)
    ProjdU, footprint = reproject_interp(TempFit[5], RefFit[5].header) # Stokes dU
    ProjdU_fit = fits.ImageHDU(data=ProjdU, header=RefFit[5].header)
    ProjMask, footprint = reproject_interp(TempFit[6], RefFit[6].header) # Stokes dU
    ProjMask_fit = fits.ImageHDU(data=ProjMask, header=RefFit[6].header)

    # Appending the Stokes I, Q, U data to the new fits file
    ProjFit = fits.HDUList()
    ProjFit.append(ProjI_fit)
    ProjFit.append(ProjdI_fit)
    ProjFit.append(ProjQ_fit)
    ProjFit.append(ProjdQ_fit)
    ProjFit.append(ProjU_fit)
    ProjFit.append(ProjdU_fit)
    # Appending the image mask
    ProjFit.append(ProjMask_fit)

    # Calculating the polarization properties
    StokesI = np.copy(ProjI)
    dI = np.copy(ProjdI)
    StokesQ = np.copy(ProjQ)
    dQ = np.copy(ProjdQ)
    StokesU = np.copy(ProjU)
    dU = np.copy(ProjdU)
    # Biased polarized intensity 
    PI_biased = (StokesQ**2.0 + StokesU**2.0)**0.5
    # PI uncertainties
    dPI = ((StokesQ*dQ)**2.0 + (StokesU*dU)**2.0)**0.5/PI_biased 
    # Mapping where the intensity is smaller than the uncertainty
    pimask = np.where(PI_biased < dPI)
    dPI_debias = np.copy(dPI)
    dPI_debias[pimask] = 0.0 #np.nan
    # De-biased polarized intensity PI 	
    PolPI = (PI_biased**2.0 - dPI_debias**2.0)**0.5  
    PolPI[pimask] = 0.0

    # Polarization fraction P
    P_biased = 100.0*PI_biased/StokesI
    PolP = 100.0*PolPI/StokesI
    # Uncertainty of polarization fraction
    dP = np.absolute(P_biased*((dPI/PI_biased)**2.0 + (dI/StokesI)**2.0)**0.5)
    # Eliminating unphysical results (setting values to 0.0 or 100.0)
    P_biased[np.where(P_biased < 0.0)] = 0.0
    P_biased[np.where(P_biased > 100.0)] = 100.0
    PolP[np.where(PolP < 0.0)] = 0.0
    PolP[np.where(PolP > 100.0)] = 100.0
    dP[np.where(dP > 100.0)] = 100.0

    # Polarization angle O
    PolO = (0.5*180.0/math.pi)*np.arctan2(StokesU, StokesQ)
    # Uncertainties dO
    dO = (0.5*180.0/math.pi)*((StokesQ*dU)**2.0 + 
                                (StokesU*dQ)**2.0)**0.5/PI_biased**2.0
    # Polarization angle B
    PolB = PolO + 90.0
    PolB[np.where(PolB > 90.0)] = PolB[np.where(PolB > 90.0)] - 180.0

    # Adding the remaining extensions
    # Polarization fraction
    ProjFit.append(fits.ImageHDU(data=P_biased, header=RefFit[7].header))
    ProjFit.append(fits.ImageHDU(data=PolP, header=RefFit[8].header))
    ProjFit.append(fits.ImageHDU(data=dP, header=RefFit[9].header))
    # Polarization angle
    ProjFit.append(fits.ImageHDU(data=PolO, header=RefFit[10].header))
    ProjFit.append(fits.ImageHDU(data=PolB, header=RefFit[11].header))
    ProjFit.append(fits.ImageHDU(data=dO, header=RefFit[12].header))
    # Polarized intensity
    ProjFit.append(fits.ImageHDU(data=PI_biased, header=RefFit[13].header))
    ProjFit.append(fits.ImageHDU(data=dPI, header=RefFit[14].header))
    ProjFit.append(fits.ImageHDU(data=PolPI, header=RefFit[15].header))
    # Adding catalog at the end
    ProjFit.append(TempFit[16])
    ProjFit.append(TempFit[17])
    # Returning the finalize HDUlist
    return ProjFit

# Function to round error decimals up
#  See https://kodify.net/python/math/round-decimals/
def round_decimals_up(number:float, decimals:int=2):
    """
    Returns a value rounded up to a specific number of decimal places.
    """
    if not isinstance(decimals, int):
        raise TypeError("decimal places must be an integer")
    elif decimals < 0:
        raise ValueError("decimal places has to be 0 or more")
    elif decimals == 0:
        return math.ceil(number)

    factor = 10 ** decimals
    return math.ceil(number * factor) / factor